% Kerr heat engine? looking for optimal parameters

clear
global nbar

nbar=2;
chimin=0.02;
chimax=2.5;
dchi=0.01;
chis=chimin:dchi:chimax;
alphamin=0.1;
alphamax=0.3;
dalpha=0.005;
alphas=alphamin:dalpha:alphamax;

velik=20;

[Alpha,Chi]=meshgrid(alphas,chis);

G=cos(Alpha).^2.*sin(Alpha).^2.*sin(Chi)/nbar.^2./((1/nbar+2*sin(Alpha).^2.*sin(Chi/2).^2).^2+sin(Alpha).^4.*sin(Chi).^2).^2;
nplus=cos(Alpha).^2*nbar+G;

Gp=cos(Alpha).^2.*sin(Alpha).^2.*Chi/nbar.^2./((1/nbar+2*sin(Alpha).^2.*Chi.^2/4).^2+sin(Alpha).^4.*Chi.^2).^2;
nplusp=cos(Alpha).^2*nbar+Gp;

%% koherentni stavy
%G=nbar*cos(Alpha).^2.*exp(-4*sin(Alpha).^2*nbar.*sin(Chi/2).^2).*besselj(1,2*sin(Alpha).^2*nbar.*sin(Chi));
%nplus=cos(Alpha).^2*nbar+G;

figure(1)
contour(alphas,chis,nplus)

figure(2)
mesh(nplusp)



function y=Dolik(x)
 global nbar
 alpha=x(1);
 chi=x(2);
 G=cos(alpha)^2*sin(alpha)^2*sin(chi)/nbar^2/((1/nbar+2*sin(alpha)^2*sin(chi/2)^2)^2+sin(alpha)^4*sin(chi)^2)^2;
 y=-cos(alpha)^2*nbar-G;
endfunction

function y=DolikCoh(x)
 global nbar
 alpha=x(1);
 chi=x(2);
 % koherentni stavy
 G=nbar*cos(alpha).^2.*exp(-4*sin(alpha).^2*nbar.*sin(chi/2).^2).*besselj(1,2*sin(alpha).^2*nbar.*sin(chi));
 y=-cos(alpha)^2*nbar-G;
endfunction



nmin=1.1;
%nmax=20;
%dn=0.1;
nmax=2000;
dn=10;
ns=nmin:dn:nmax;
nsiz=size(ns,2);
alphopt=zeros(1,nsiz);
chiopt=zeros(1,nsiz);
alphoptC=zeros(1,nsiz);
chioptC=zeros(1,nsiz);
energopt=zeros(1,nsiz);
energoptCoh=zeros(1,nsiz);
arpg0=[0.25;1.2];
arpg=arpg0;
arpgC=arpg;
tic
for n=1:nsiz
 nbar=ns(n);
 vals=fminunc('Dolik',arpg);
 alphopt(n)=vals(1);
 chiopt(n)=vals(2);
 energopt(n)=-Dolik(vals);
 arpg=vals*(1+0.01*randn);
 % Coherent
 arpgC=[vals(1)*1.26;vals(2)];
 vals=fminunc('DolikCoh',arpgC);
 alphoptC(n)=vals(1);
 chioptC(n)=vals(2);
 energoptCoh(n)=-DolikCoh(vals);
% arpgC=vals*(1+0.01*randn);
end
toc

enrelat=energopt./ns;
enrelatCoh=energoptCoh./ns;

enrelMaxCoh=1.5819;
enrelMaxTh=1.325;

figure(3)
plot(ns,[alphopt;alphoptC])
set(gca,'FontSize',1.5*velik)


figure(4)
plot(ns,[chiopt;chioptC])
set(gca,'FontSize',1.5*velik)


figure(5)
plot(ns,[energopt;energoptCoh])
set(gca,'FontSize',velik)


figure(6)
plot(ns,enrelat)
hold on
%plot(ns,enrelatCoh,'r')
plot([nmin,nmax],enrelMaxTh*[1,1],'--')
%plot([nmin,nmax],enrelMaxCoh*[1,1],'r--')
hold off
set(gca,'FontSize',2*velik)
%set(gca,'FontSize',velik)

figure(7)
plot(ns,alphoptC-alphopt)

uncn=zeros(1,nsiz);
uncTerm=zeros(1,nsiz);
enCoh=zeros(1,nsiz);
auxP=zeros(4,nsiz);
for n=1:nsiz
 nbar=ns(n);
 s=sin(alphopt(n));
 c=cos(alphopt(n));
 chi=chiopt(n);
 nmean=energopt(n);
 zavorka1=1/nbar+2*s^2*sin(chi/2)^2;
 ssinch1=s^2*sin(chi);
 pomA=4*ssinch1*zavorka1/nbar^2/(zavorka1^2+ssinch1^2)^3;
 zavorka2=1/nbar+2*s^2*sin(chi)^2;
 ssinch2=s^4*sin(2*chi)^2;
 pomB=ssinch2/nbar^2/(zavorka2^2+ssinch2)^3;
 varian=c^4*(2*nbar^2+pomA+pomB)+nmean-nmean^2;
 uncn(n)=sqrt(varian);
 uncTerm(n)=sqrt(nmean^2+nmean);
 auxP(1,n)=2*nbar^2;
 auxP(2,n)=pomA;
 auxP(3,n)=pomB;
 auxP(4,n)=nmean;
% uncn(n)=pomA;
% Coherent states
% term1=2*exp(-4*s^2*nbar*sin(chi/2)^2)*besselj(1,2*nbar*s^2*sin(chi));
% term2=1/2*exp(-4*s^2*nbar*sin(chi)^2)*besselj(2,2*nbar*s^2*sin(2*chi));
% adagadagaa=c^4*nbar^2*(3/2+term1+term2);
% varian=adagadagaa+nmean-nmean^2;
% uncn(n)=sqrt(varian);
% uncTerm(n)=sqrt(nbar^2/2+nbar);  % for coherent states with random phases
% enCoh(n)=nbar;
end

figure(7)
plot(ns,uncn)
plot(ns,energopt)
hold on
plot(ns,energopt-uncn,'--')
plot(ns,energopt+uncn,'--')
%plot(ns,enCoh,'r')
%plot(ns,enCoh-uncTerm,'r--')
%plot(ns,enCoh+uncTerm,'r--')
hold off
set(gca,'FontSize',velik)

figure(8)
%plot(ns,auxP)
%plot(ns,uncn./uncTerm)
%plot(ns,energopt./enCoh)
plot(ns,uncn)
set(gca,'FontSize',velik*1.5)

chi=0.84
alpha=0.244
nbar=5
G=cos(alpha).^2.*sin(alpha).^2.*sin(chi)/nbar.^2./((1/nbar+2*sin(alpha).^2.*sin(chi/2).^2).^2+sin(alpha).^4.*sin(chi).^2).^2;
nplusS=cos(alpha).^2*nbar+G
nminusS=cos(alpha).^2*nbar-G
